Análisis de Regresión con R

Comparación de modelos de regresión con R

Artículo científico

Datos

• Datos tomados del artículo referenciado previamente.
• Fuente de datos - Artículo en PLOS | ONE.
• En la base de datos se tienen dos grupos bajo análisis:
  • Grupo Long jumpers: id 1 a 23.
  • Grupo untrained men: id 24-45.

Lectura de datos

library(readxl)
library(tidyverse)
datos <- read_xlsx("../data/Data_Paper_Plos_One_Muscle.xlsx", skip = 3,
                   na = "N/A", n_max = 47) %>% 
  rename(RA_takeof_leg = RA...7,
         OB_takeof_leg = OB...8,
         PM_takeof_leg = PM...9,
         QL_takeof_leg = QL...10,
         ES_takeof_leg = ES...11,
         Gmax_takeof_leg = Gmax...12,
         Gmed_takeof_leg = Gmed...13,
         IL_takeof_leg = IL...14,
         RA_free_leg = RA...15,
         OB_free_leg = OB...16,
         PM_free_leg = PM...17,
         QL_free_leg = QL...18,
         ES_free_leg = ES...19,
         Gmax_free_leg = Gmax...20,
         Gmed_free_leg = Gmed...21,
         IL_free_leg = IL...22,
         id = ID,
         edad = `Age (years)`,
         altura_cm = `Height (cm)`,
         imc = `Body mass (kg)`,
         dist_salto_cm = `long jump distance (cm)`,
         sprint_100m_seconds = `100-m sprint time (s)`,
         grasa_subcut_cm2 = `Subcutaneous fat CSA (absolute value, cm2)`) %>% 
  mutate(type = if_else(id %in% c(1:23), true = "Long jumpers",
                        false = "Untrained men"))
datos

Objetivos

• Replicar análisis estadísticos aplicados en el artículo científico de interés.
• Evidenciar la relación existente entre características anatómicas de atletas vs rendimiento en salto largo.
• Evaluar otros métodos de statistical learning y compararlos con los resultados obtenidos por los autores.

Resultados del paper

Correlaciones

• Aunque fueron numerosos los resultados obtenidos por los autores, para el objetivo de este documento se destacan los siguientes:
  • La relación entre el área transversal relativa (CSA) del recto abdominal (AR) del lado de la pierna de despegue y el mejor registro personal para el salto largo.
    • Correlación: 0.674
    • Valor p: 0.004 (estadísticamente significativo)
• Las correlaciones (con intervalo de confianza del 95%) se presentan en la siguiente tabla:

Gráfico de dispersión

Original

Réplica con R

library(ggplot2)
datos %>% 
  ggplot(data = ., aes(x = RA_takeof_leg, y = dist_salto_cm)) +
  geom_point(size = 3) +
  labs(x = expression('Relative CSA of RA takeoff leg side - cm'^"2"/'kg'^"2/3"),
       y = "Personal best record of long jump (cm)") +
  geom_smooth(method = "lm", se = FALSE, lty = 3, lwd = 1, color = "black") +
  theme_light()

Predichos vs Reales (paper)

Resultados adicionales con R

Distribuciones

datos %>% 
  select_if(is.numeric) %>% 
  select(-id) %>% 
  gather(key = "variable", value = "valor") %>% 
  ggplot(data = ., aes(x = valor)) +
  facet_wrap(facets = ~variable, scales = "free", ncol = 4) +
  geom_histogram(aes(y = ..density..), bins = 10, color = "black", 
                 fill = "gray60") +
  geom_density(fill = "gray50", alpha = 0.18) +
  geom_rug() +
  labs(x = "", y = "Densidad") +
  theme_light() +
  theme(strip.background = element_rect(fill = "deepskyblue4"),
        strip.text = element_text(color = "black"))

Gráficos cuantil cuantil

library(qqplotr)
datos %>% 
  select_if(is.numeric) %>% 
  select(-id) %>% 
  gather(key = "variable", value = "valor") %>% 
  ggplot(data = ., aes(sample = valor)) +
  facet_wrap(facets = ~variable, scales = "free", ncol = 4) +
  geom_qq_band(fill = "gray25") +
  stat_qq_line(color = "darkgreen") +
  stat_qq_point(color = "black", size = 0.8) +
  labs(x = "Cuantiles teóricos", y = "Cuantiles muestrales") +
  theme_light() +
  theme(strip.background = element_rect(fill = "deepskyblue4"),
        strip.text = element_text(color = "black"))

Comparativos

• Se comparan registros de cross-sectional area (CSA) de la pierna de despeque (takeoff) vs la pierna libre (free). Las variables (músculos) a comparar son las siguientes:
  • RA: recto abdominal.
  • OB: oblicuos internos y externos.
  • PM: psoas mayor.
  • QL: cuadrado lumbar.
  • ES: erector spinae.
  • Gmax: gluteo mayor.
  • Gmed: gluteos medio y mínimo.
  • IL: iliaco

df_takeoff_leg <- datos %>% 
  select(RA_takeof_leg:IL_takeof_leg) %>% 
  gather(key = "variable", value = "valor") %>% 
  mutate(tipo = "TakeoffLeg")

df_free_leg <- datos %>% 
  select(RA_free_leg:IL_free_leg) %>% 
  gather(key = "variable", value = "valor") %>% 
  mutate(tipo = "FreeLeg")

df_takeoff_free <- df_takeoff_leg %>% 
  bind_rows(df_free_leg)

df_takeoff_free %>% 
  separate(col = variable, into = c("variable", "v1", "v2")) %>% 
  select(-c(v1, v2))  %>% 
  ggplot(data = ., aes(x = tipo, y = valor, fill = tipo)) +
  facet_wrap(facets = ~variable, scales = "free", ncol = 4) +
  geom_boxplot(color = "black") +
  scale_fill_manual(values = c("darkgreen", "gold4")) +
  labs(x = "Tipo de pierna", y = "") +
  theme_light() +
  theme(strip.background = element_rect(fill = "deepskyblue4"),
        strip.text = element_text(color = "black"),
        legend.position = "none")

Shapiro Wilk

Se comprueba la normalidad de las variables (\(\alpha = 0.05\)), bajo el siguiente juego de hipótesis:

\[H_0: X \sim N(\mu, \sigma^2)\\ H1: x \nsim N(\mu, \sigma^2)\]

datos %>% 
  select_if(is.numeric) %>% 
  select(-id) %>% 
  gather(key = "variable", value = "valor") %>% 
  group_by(variable) %>% 
  summarise(valor = list(valor)) %>% 
  ungroup() %>% 
  group_by(variable) %>% 
  mutate(shapiro_valorP = round(shapiro.test(unlist(valor))$p.value, digits = 5),
         Resultado = if_else(shapiro_valorP <= 0.05, true = "No normalidad",
                             false = "Sí normalidad"))  %>% 
  select(-valor)

• Nota: aunque los autores mencionan que fueron aplicadas las pruebas de Shapiro Wilk para comprobar el supuesto de normalidad, en la tabla anterior se evidencia que algunas variables (ej. distancia de salto, grasa subcutanea, entre otras) no se distribuyen de forma normal. Este resultado tiene connotaciones de importancia, ya que las correlaciones podrían ser obtenidas mediante métodos no paramétricos.

Matriz de correlaciones

• Se construye la matriz de correlaciones (método de Pearson).
• La variable que presente mayor correlación lineal con la longitud del salto, será tenida en cuenta para estructurar el modelo de regresión lineal simple (RLS). Con las demás variables se construye el modelo de regresión lineal múltiple (RLM).

library(corrplot)
library(RColorBrewer)
datos %>% 
  select_if(is.numeric) %>% 
  select(-id) %>% 
  cor(use = "complete.obs") %>% 
  corrplot(method = "pie",
           type = "upper",
           diag = FALSE,
           tl.cex = 0.8,
           tl.srt = 45,
           addgrid.col = "black",
           order = "hclust",
           col = brewer.pal(n = 10, name = "Spectral"))

• Notas:
  • Se evidencia alta correlación entre algunas variables. Este patrón sugiere problemas de multicolinealidad al ajustar un modelo de RLM.
  • Atletas con mayor altura tienden a presentar saltos de menor distancia. Este comportamiento es obtenido de igual forma en aquellos que aumentan la velocidad del sprint al saltar, es decir, que a mayor velocidad, menor longitud de salto.
  • Atletas con el músculo recto abdominal de mayor longitud, presentan mayores distancias en sus saltos. Con esta variable se ajusta el modelo de RLS.

Regresión Lineal Simple (RLS)

• Se propone cuantificar la unidad de cambio en la distancia del salto vs la longitud (cm²/kg^2/3) del músculo recto abdominal (RA), para la pierna de despegue (takeoff).
• El modelo de RLS puede ser expresado de la siguiente manera: \[Y = \beta_0\ + \beta_1X \\ \hat{Y} = \hat{\beta_0}\ + \hat{\beta_1}X\ + \epsilon\]
• Escrito de cara al fenómeno bajo estudio, el modelo queda expresado como sigue: \[Distancia = \hat{\beta_0}\ + (\hat{\beta_1}\times RA) + \epsilon\]
• Este modelo de RLS es ajustado a través del Método de Mínimos Cuadrados.

Modelo Lineal con lm()

mod_rls <- lm(dist_salto_cm ~ RA_takeof_leg, data = datos)
summary(mod_rls)

Call:
lm(formula = dist_salto_cm ~ RA_takeof_leg, data = datos)

Residuals:
   Min     1Q Median     3Q    Max 
-72.26 -22.41   8.38  22.24  65.10 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    589.584     49.894  11.817 9.66e-11 ***
RA_takeof_leg   11.928      4.443   2.685   0.0139 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 40.89 on 21 degrees of freedom
  (22 observations deleted due to missingness)
Multiple R-squared:  0.2555,    Adjusted R-squared:  0.2201 
F-statistic: 7.208 on 1 and 21 DF,  p-value: 0.01387

• Nota: el resultado anterior sugiere que la variable RA_takeof_leg es estadísticamente significativa (\(valor\ p =0.0139\)) sobre la variabilidad observada en la distancia del salto. Además, se puede inferir que por cada unidad que aumenta RA_takeof_leg, la distancia de salto es 11.928 centímetros mayor. La variable RA_takeof_leg explica 22.01% de la variabilidad observada en la distancia de salto.

Significancia Estadística

anova(mod_rls)

Analysis of Variance Table

Response: dist_salto_cm
              Df Sum Sq Mean Sq F value  Pr(>F)  
RA_takeof_leg  1  12050 12050.0  7.2078 0.01387 *
Residuals     21  35108  1671.8                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residuales

par(mfrow = c(2, 2))
plot(mod_rls)

• Normalidad (Shapiro Wilk) de los residuales:

shapiro.test(residuals(mod_rls))

    Shapiro-Wilk normality test

data:  residuals(mod_rls)
W = 0.96191, p-value = 0.5029

• Heterocedasticidad - Breusch Pagan Test:

library(lmtest)
bptest(mod_rls)

    studentized Breusch-Pagan test

data:  mod_rls
BP = 0.0025838, df = 1, p-value = 0.9595

• Autocorrelación - Test de Durbi Watson:

dwtest(mod_rls)

    Durbin-Watson test

data:  mod_rls
DW = 2.1677, p-value = 0.7123
alternative hypothesis: true autocorrelation is greater than 0

• Se comprueba que existe normalidad de los residuos, son homocedasticos y no existe autocorrelación de los mismos.

Reales vs Predichos

reales <- datos$dist_salto_cm[!is.na(datos$dist_salto_cm)]
predichos_rls <- mod_rls$fitted.values
data.frame(
  Real = reales,
  Predichos = predichos_rls
) %>% 
  ggplot(data = ., aes(x = Predichos, y = Real)) +
  geom_point() +
  labs(x = ("Valores predichos de distancia (cm)"),
       y = "Valores reales de distancia (cm)") +
  geom_smooth(method = "lm", se = FALSE, lty = 3, lwd = 1, color = "black") +
  theme_light()

• La correlación entre los valores reales y valores predichos por el modelo de RLS es: 0.505; más baja que la reportada por los autores al ajustar un modelo de RLM, cuyo valor es igual a 0.892

Regresión Lineal Múltiple (RLM)

• Para la construcción del modelo de RLM se comprueba la multicolinealidad de las variables y se proponen cuatro alternativas:
  • Modelo 0: modelo con todas las variables incluidas como predictoras.
  • Modelo 1: Modelo de RLM con eliminación de variables por valores de correlación y factor inflacionario de varianza - VIF..
  • Modelo 2: Modelo de RLM con eliminación de variables por VIF a través del métodos Stepwise (utilizado por los autores del artículo).
  • Modelo 3: Modelo de regresión por componentes principales (mínimos cuadrados parciales).
• Los modelos son comparados a través del R² ajustado (mayor mejor) y el cuadrado medio del error (CME).

Modelo 0

mod_rlm0 <- lm(dist_salto_cm ~ .,
               data = datos %>% select_if(is.numeric) %>% select(-id))
summary(mod_rlm0)

Call:
lm(formula = dist_salto_cm ~ ., data = datos %>% select_if(is.numeric) %>% 
    select(-id))

Residuals:
       1        2        3        4        5        6        7        8        9 
-0.24001 -1.43196 -0.91660 -0.75285 -0.04370  1.48290 -0.44264  1.69525 -0.39388 
      10       11       12       13       14       15       16       17       18 
 1.11512 -1.41008  0.27780 -0.83303 -0.28836  1.11155 -0.39267  0.08616  1.18270 
      19       20       21       22       23 
 0.62738 -0.80017  1.33190 -1.93656  0.97174 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)  
(Intercept)         1266.6968   182.3241   6.948   0.0910 .
edad                 -23.4171     3.4665  -6.755   0.0936 .
altura_cm              6.1246     1.3465   4.549   0.1378  
imc                   -2.7528     1.3980  -1.969   0.2992  
sprint_100m_seconds  -92.8684     5.1288 -18.107   0.0351 *
RA_takeof_leg        124.1035    12.0047  10.338   0.0614 .
OB_takeof_leg          7.1425     2.1907   3.260   0.1895  
PM_takeof_leg          9.0525     2.7581   3.282   0.1883  
QL_takeof_leg          7.8868     1.4157   5.571   0.1131  
ES_takeof_leg         35.5716     5.3988   6.589   0.0959 .
Gmax_takeof_leg       -1.0174     0.8984  -1.132   0.4605  
Gmed_takeof_leg      -10.5371     0.9414 -11.193   0.0567 .
IL_takeof_leg         -0.4193     2.4356  -0.172   0.8915  
RA_free_leg          -86.7869     9.1217  -9.514   0.0667 .
OB_free_leg           -2.4849     1.2159  -2.044   0.2897  
PM_free_leg          -17.9138     3.3439  -5.357   0.1175  
QL_free_leg           -7.7265     3.9823  -1.940   0.3030  
ES_free_leg          -49.2399     6.2703  -7.853   0.0806 .
Gmax_free_leg          2.4543     0.5668   4.330   0.1445  
Gmed_free_leg         12.7940     0.9606  13.319   0.0477 *
IL_free_leg            2.8488     2.3869   1.194   0.4440  
grasa_subcut_cm2      -2.3272     0.3648  -6.379   0.0990 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.83 on 1 degrees of freedom
  (22 observations deleted due to missingness)
Multiple R-squared:  0.9995,    Adjusted R-squared:  0.9891 
F-statistic: 96.23 on 21 and 1 DF,  p-value: 0.08023

• El R² ajustado es próximo a 1, sin embargo, pocas variables independientes (predictoras) son estadísticamente significativas (\(p<0.05\)), patron de comportamiento que sugiere problemas de multicolinealidad. La siguiente gráfica muestra anomalías en los residuales del modelo.

par(mfrow = c(2, 2))
plot(mod_rlm0)

• Estos resultados evidencian que el modelo ajustado no es adecuado para explicar la variabilidad observada en la distancia de salto. Es importante resaltar que el R² no es la mejor medida de bondad de ajuste, ya que éste aumenta en función del número de predictores.

# Valores predichos por mod_rlm0
predichos_rlm0 <- mod_rlm0$fitted.values 

Modelo 1

• Selección de predictores teniendo en cuenta valores iguales o inferiores a 0.70 de correlación. La selección es posible hacerla con la bibliteca usdm que posee funciones flexibles para diagnósticos de multicolinealidad.

# Variables predictoras
df_predictoras <-datos %>%
  select_if(is.numeric) %>% dplyr::select(-c(id, dist_salto_cm)) %>% 
  as.matrix()

library(usdm)
vifcor(x = df_predictoras, th = 0.70)

10 variables from the 21 input variables have collinearity problem: 
 
ES_free_leg RA_takeof_leg OB_takeof_leg PM_takeof_leg Gmed_takeof_leg imc Gmax_free_leg Gmax_takeof_leg ES_takeof_leg QL_free_leg 

After excluding the collinear variables, the linear correlation coefficients ranges between: 
min correlation ( OB_free_leg ~ sprint_100m_seconds ):  -0.001734317 
max correlation ( IL_free_leg ~ IL_takeof_leg ):  0.6817036 

---------- VIFs of the remained variables -------- 

• Este resultado muestra que de las 21 variables consideradas en el modelo inicial (modelo 0) 10 de ellas presentan problemas de colinealidad. Tomando un límite de 0.70 de correlación como criterio de exclusión de predictores, la correlación máxima presente en las variables seleccionadas, es de 0.68. La tabla anterior muestra los predictores que podrían hacer parte del modelo.

• Nuevo modelo:

mod_rlm1 <- lm(dist_salto_cm ~ edad + altura_cm + sprint_100m_seconds
               + QL_takeof_leg + IL_takeof_leg + RA_free_leg + OB_free_leg
               + PM_free_leg + Gmed_free_leg + IL_free_leg + grasa_subcut_cm2,
               data = datos)
summary(mod_rlm1)

Call:
lm(formula = dist_salto_cm ~ edad + altura_cm + sprint_100m_seconds + 
    QL_takeof_leg + IL_takeof_leg + RA_free_leg + OB_free_leg + 
    PM_free_leg + Gmed_free_leg + IL_free_leg + grasa_subcut_cm2, 
    data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.024 -13.760  -0.812  13.552  38.971 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         1859.2274   339.7523   5.472 0.000194 ***
edad                   4.8879     5.8586   0.834 0.421854    
altura_cm             -3.8708     1.5704  -2.465 0.031409 *  
sprint_100m_seconds  -58.9174    20.8987  -2.819 0.016692 *  
QL_takeof_leg         -2.5243     5.5584  -0.454 0.658553    
IL_takeof_leg         10.7345     6.0471   1.775 0.103514    
RA_free_leg            6.0558     5.5662   1.088 0.299875    
OB_free_leg           -0.4214     2.0322  -0.207 0.839505    
PM_free_leg           -0.1120     3.1732  -0.035 0.972481    
Gmed_free_leg         -0.6212     1.5826  -0.393 0.702176    
IL_free_leg           -4.6031     6.7318  -0.684 0.508261    
grasa_subcut_cm2       0.4789     0.8361   0.573 0.578311    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 29.51 on 11 degrees of freedom
  (22 observations deleted due to missingness)
Multiple R-squared:  0.7968,    Adjusted R-squared:  0.5936 
F-statistic: 3.921 on 11 and 11 DF,  p-value: 0.01622

• El R² ajustado es igual a 0.5936 y sólo son estadísticamente significativas las variables altura y sprint. En general, el modelo es estadísticamente significativo (\(p = 0.01622\)) para explicar la variabilidad en la distancia de salto.

# Valores predichos por mod_rlm1
predichos_rlm1 <- mod_rlm1$fitted.values

Modelo 2

Modelo 3

Modelo 4